function output = LikelihoodFunction_Het_v14(parameters, TreatmentPop, RC, J, Popsize,f,d,mFixed,mFlexMin,Sim,acceptsContractData,commitmentAmountData,attendsCampData,hbp_indData,sickness_indData, employeeData, taking_medsData, genderData, literateData)

%%%%%%%%%%%%%% Set parameters from input matrix in optimizer %%%%%
 
betaMean	=	parameters(1);  	 	% Mean of present bias parameter, 0<beta<=1
betaSD = parameters(2);  				% SD of present bias parameter, 0<beta<=1
betaHatBeta1	=	parameters(3);   	% Parameter 1 for Beta Distribution, defining distance between beta and 1 for betaHat (level of sophistication), beta <=betaHat
betaHatBeta2	=	parameters(4);   	% Parameter 2 for Beta Distribution, defining distance between beta and 1 for betaHat (level of sophistication), beta <= betaHat
delta = 	parameters(5);      		% Exponential discount factor (only applies to delayed health benefit). This is important! No exponential discounting from t=0 to t=1, since only 1 week. Can be changed.
cMean	=	parameters(6);      		% Utility cost Mean of attending health camp
cSD = parameters(7);					% Utility cost SD of attending health camp (in basic specification, individual knows draw at time of contract choice)
eps1SD = parameters(8);					% Epsilon draw from mean 0 distribution, SD estimated, for t=0 perceptions of t=1 attendence value 
base_b	=	parameters(9);          		% Delayed health benefit of attending [Bringing in as data right now, but cuold bring into estimation later]
bObsHet_hbp= parameters(10);				% Shift in mean of benefit heterogeneity by indicator for severity of hbp
bObsHet_sk = parameters(11);				% Shift in mean of benefit heterogeneity by indicator for sickness
bObsHet_sx = parameters(12);				% Shift in mean of benefit heterogeneity by indicator for gender
btObsHet_sx = parameters(13);				% Shift in mean of benefit heterogeneity by indicator for gender
bthObsHet_sx = parameters(14);				% Shift in mean of benefit heterogeneity by indicator for gender
cObsHet_emp= parameters(15);				% Shift in mean of cost heterogeneity by whether resp is laborer or employee
bObsHet_med= parameters(16);				% Shift in mean of benefit heterogeneity by whether resp takes bp meds
bObsHet_lit= parameters(17);				% Shift in mean of benefit heterogeneity by whether resp takes bp meds


%%% Translate into individual-level preferences for population, using random coefficients approach, Sim simulations per individual 

betaPopTemp = betaMean*ones(1,Sim); 
betaPopSD = (RC(:,1)*betaSD)';
betaPop = ones(Popsize,1)*(betaPopTemp + betaPopSD);
betaPop =  betaPop+(btObsHet_sx*genderData);
betaPop = min(betaPop,1);
betaPop = max(betaPop,0);

InvBeta = betainv(J, betaHatBeta1,betaHatBeta2);
betaHatPopTemp = ones(Popsize,1)*InvBeta;
betaHatPopTemp = betaHatPopTemp + (genderData*bthObsHet_sx);
betaHatPop =  betaPop + (1-betaPop).* betaHatPopTemp;

costPopTemp = cMean*ones(1,Sim);
costPopSD = (RC(:,2)*cSD)';
costPop = (ones(Popsize,1)*(costPopTemp + costPopSD))  + (cObsHet_emp*employeeData);
costPop = max(costPop,0);

eps1Pop = ones(Popsize,1)*(RC(:,3)*eps1SD)';
eps2Pop = eps1Pop ;

%%%benefit b with shifters
btemp = base_b*ones(1,Sim);
b = ((ones(Popsize,1)*btemp)+(bObsHet_hbp*hbp_indData)+(bObsHet_sk*sickness_indData)+(bObsHet_med*taking_medsData)+(bObsHet_sx*genderData)+(bObsHet_lit*literateData));

%%% Compute key conditions from theory model to help simplify choice computation in latter part of likelihood function 

% First, compute whether t=0 (long run) self wants to attend to the camp in the future

wantsToAttend = (costPop+f <= delta*(b + eps1Pop)); 
wantsToAttendWithDiscount = (costPop+f-d <= delta*(b + eps1Pop)); 

% Now, compute whether t=0 self believes she wlll attend at t=1, without any commitment contract or discount 

believesWillAttendWithoutCommit =  (costPop+f <= betaHatPop*delta.*(b + eps1Pop));  

% Now compute whether t=0 self has positive demand for *some* possible commitment contract (maybe not the one offered)

tasteForCommitment = zeros(Popsize,Sim);

for i = 1:Popsize
  for j = 1:Sim 
  
	if (wantsToAttend(i,j)==1 && believesWillAttendWithoutCommit(i,j)==0)
		tasteForCommitment(i,j)=1;
	else
		tasteForCommitment(i,j)=0;    
	end
	
  end
end

% Now compute perceived optimal commitment amount for each individual. This is just the minimum amount that ensures they will attend the camp ex post, given 
% their beliefs about their beta, betahat

for i = 1:Popsize
	for j = 1:Sim
		mStar(i,j) = costPop(i,j) - betaHatPop(i,j)*delta*(b(i,j)+eps1Pop(i,j));
		mFixLargeEnough(i,j) = (mFixed > mStar(i,j));	
	end	
end 

acceptsContract = zeros(Popsize,Sim);
acceptsContract(:,:)=NaN;
acceptsContractLikelihood = zeros(Popsize,Sim);
commitmentAmountChosen = zeros(Popsize,Sim);
commitmentAmountChosen(:,:)=NaN;
commitmentAmountChosenLikelihood = zeros(Popsize,Sim);

% =============================================================
% Treatments, Contract Choices and Commitment Amount Chosen
% =============================================================

for i = 1:Popsize
		for j = 1:Sim
	
	% In our model right now, if taste for commitment is 0, we assume they are indifferent and don't choose commitment contract. Simplification, will need to fix later.  	
	
	if TreatmentPop(i,1) == 2  % FIXED COMMITMENT
		if (tasteForCommitment(i,j)==1 && mFixLargeEnough(i,j)==1)  % Taste for commitment + fixed commitment strong enough
			acceptsContract(i,j)= 1 - max(min((0.1/(costPop(i,j) + f - betaHatPop(i,j)*delta*(b(i,j)+eps1Pop(i,j)))),0.01),0); % Second term is smoothing function related to "BelievesWillAttendWitoutCommit"
			commitmentAmountChosen(i,j) = mFixed;
		else
			acceptsContract(i,j)= 0 - max(min((-0.1/(costPop(i,j) + f - betaHatPop(i,j)*delta*(b(i,j)+eps1Pop(i,j)))),0),-0.01); % Second term is smoothing related to wants to Attend, likely = 0.  
		end																												 	
		
		
	elseif TreatmentPop(i,1) == 3  % FIXED COMMITMENT WITH DISCOUNT
		if (wantsToAttendWithDiscount(i,j)==1 && mFixLargeEnough(i,j)==1)  % Wants to go + predicts will go with commitment / discount (note no need to have taste for commitment)
			acceptsContract(i,j)= 1  - max(min((0.1/(costPop(i,j) + f - d - delta*(b(i,j)+eps1Pop(i,j)))),0.01),0); %%% Smoothing, related to Wants To Attend with Discount Condition
			commitmentAmountChosen(i,j) = mFixed;
		else   % Either doesn't want to go at all, or predicts would not go even with commitment amount
			acceptsContract(i,j)=0  - max(min((-0.1/(costPop(i,j) + f - d - delta*(b(i,j)+eps1Pop(i,j)))),0),-0.01); %%% Smoothing, related to Wants To Attend with Discount Condition
		end   
		
	elseif TreatmentPop(i,1) == 4  % FLEXIBLE COMMITMENT
		if (tasteForCommitment(i,j)==1)  % If contract is flexible, any taste for commitment will suffice to sign up
			acceptsContract(i,j)= 1 - max(min((0.1/(costPop(i,j) + f - betaHatPop(i,j)*delta*(b(i,j)+eps1Pop(i,j)))),0.01),0); % Second term is smoothing function related to "BelievesWillAttendWitoutCommit";
			commitmentAmountChosen(i,j) = max(mStar(i,j), mFlexMin);
		else   % No reason to sign up if no taste for commitment
			acceptsContract(i,j)=0 - max(min((-0.1/(costPop(i,j) + f - betaHatPop(i,j)*delta*(b(i,j)+eps1Pop(i,j)))),0),-0.01);
		end
    
	elseif TreatmentPop(i,1) == 5  % FLEXIBLE COMMITMENT WITH DISCOUNT
		if (wantsToAttendWithDiscount(i,j)==1) % Will sign up just for discount (and will pick sufficient commitment if reqd)
			acceptsContract(i,j)= 1 - max(min((0.1/(costPop(i,j) + f - d - delta*(b(i,j)+eps1Pop(i,j)))),0.01),0); %%% Smoothing, related to Wants To Attend with Discount Condition
			commitmentAmountChosen(i,j) = max(mStar(i,j), mFlexMin);
		else
			acceptsContract(i,j)=0 - max(min((-0.1/(costPop(i,j) + f - d - delta*(b(i,j)+eps1Pop(i,j)))),0),-0.01); %%% Smoothing, related to Wants To Attend with Discount Condition
		end
	end
    
	end	
end 

%%%%%%%%%%%%%%% We keep treatment 0 and 1 out of loop above, because we know for sure acceptsContract and commitment amount will be NaN, and we won't 
%%%%%%%%%%%%%%% include that in likelihood function for those people. We also keep out CommitmentAmountChosen for people who don't choose commitment contract, already filled to NaN. 

%%%%%%%%%%%%%%% Generate Attends Camp Variables 

attendsCamp = zeros(Popsize,Sim);

ControlAttend =  (costPop + f <= betaPop*delta.*(b+eps2Pop)) - min(max(0.01./(betaPop*delta.*(b+eps2Pop) - costPop - f),-0.01),0.01);  %%%%% Second term is smoothing function limiting to true value, so optimization works with discrete simulations. 
Control = (TreatmentPop==0);
ControlT = Control*ones(1,Sim);
attendsCamp = attendsCamp + ControlT.*ControlAttend;

DiscountAttend = (costPop + f - d <= betaPop*delta.*(b+eps2Pop)) - min(max(0.01./(betaPop*delta.*(b+eps2Pop) - costPop - f + d),-0.01),0.01); %%%%% Second term is smoothing function limiting to true value, so optimization works with discrete simulations.
Discount = (TreatmentPop==1);
DiscountT = Discount*ones(1,Sim);
attendsCamp = attendsCamp + DiscountT.*DiscountAttend;

CCAttend = (costPop <= betaPop*delta.*(b+eps2Pop) + commitmentAmountChosen) - min(max(0.01./(betaPop*delta.*(b+eps2Pop) + commitmentAmountChosen - costPop),-0.01),0.01); %%%%% Second term is smoothing function limiting to true value, so optimization works with discrete simulations.

T2 = (TreatmentPop==2);
T2Acc = ((TreatmentPop==2)*ones(1,Sim)).*(acceptsContractData==1);
T2Dec = ((TreatmentPop==2)*ones(1,Sim)).*(acceptsContractData==0);
T3 = (TreatmentPop==3);
T3Acc = ((TreatmentPop==3)*ones(1,Sim)).*(acceptsContractData==1);
T3Dec = ((TreatmentPop==3)*ones(1,Sim)).*(acceptsContractData==0);
T4 = (TreatmentPop==4);
T4Acc = ((TreatmentPop==4)*ones(1,Sim)).*(acceptsContractData==1);
T4Dec = ((TreatmentPop==4)*ones(1,Sim)).*(acceptsContractData==0);
T5 = (TreatmentPop==5);
T5Acc = ((TreatmentPop==5)*ones(1,Sim)).*(acceptsContractData==1);
T5Dec = ((TreatmentPop==5)*ones(1,Sim)).*(acceptsContractData==0);

CCAttend(isnan(CCAttend)) = 0;  %%%%%%%%%%%%%%%%% Replace NaN for CCAttend with 0s. This just makes computation below work, has no impact since NaNs are not counted in likelihodd below. 

attendsCamp = attendsCamp + (T2Acc + T3Acc + T4Acc + T5Acc).*CCAttend + (T2Dec + T3Dec + T4Dec + T5Dec).*ControlAttend;

%%%%%%%%% This brings in observed data, and creates likelihood portion attributed to Camp Attendence by only counting 
%%%%%%%%% Predicted Probability of Choosing Thing That Was Actually Chosen for Camp Attendence / Contract Choice. 

%%%%%%%%% Contract Choice Likelihood 

%%%%%%%%% DO NOT USE LOGICAL THAT SETS THESE THINGS EQUAL, SMOOTHING MAKES THIS NOT POSSIBLE / DESIRABLE

acceptsContractLikelihood = (sum((acceptsContractData.*acceptsContract),2) + sum((1-acceptsContractData).*(1-acceptsContract),2))/Sim;
acceptsContractLikelihood(isnan(acceptsContractLikelihood)) = 0; %%%%%%%%%%%%%%%%% Replace NaN for CCAttend with 0s. This just makes computation below work, has no impact since NaNs are not counted in likelihodd below.

%%%%%%%%% Attends Likelihood in copmutation below must be conditional on having made choice that person actually makes 

ContractCorrect1 = (acceptsContractData <= acceptsContract + 0.01);
ContractCorrect2 = (acceptsContractData >= acceptsContract - 0.01);
ContractCorrectT = ContractCorrect1 + ContractCorrect2;
ContractCorrect = (ContractCorrectT==2);
NumContractCorrect = sum(ContractCorrect,2);

attendsLikelihood = zeros(Popsize,1);
attendsLikelihood = (Control + Discount).*(sum((attendsCampData.*attendsCamp),2) + sum((1-attendsCampData).*(1-attendsCamp),2))/Sim; 
attendsLikelihood = attendsLikelihood + (T2+T3+T4+T5).*(sum((attendsCampData.*attendsCamp.*ContractCorrect),2) + sum((1-attendsCampData).*(1-attendsCamp).*ContractCorrect,2))./max(NumContractCorrect,1);

%%%%%%%%% Construct Likelihood for Commitment Amount Chosen as some density measure for actual amount chosen (within some defined window or defined kernel) 
%%%%%%%%% This is uniform kernel approach with unirofm kernel of 5 

%%%%%%%%% Commitment Amount Conditional must be conditional on actually having chosen contract in simulation 

commitmentAmtChosenLB = commitmentAmountData - 5;
commitmentAmtChosenUB = commitmentAmountData + 5;
commitmentAmtLikelihoodT1 = (commitmentAmountChosen <= commitmentAmtChosenUB); 
commitmentAmtLikelihoodT2 = (commitmentAmountChosen >= commitmentAmtChosenLB); 
commitmentAmtLikelihoodT3 = commitmentAmtLikelihoodT1 + commitmentAmtLikelihoodT2;
commitmentAmtLikelihoodT = (commitmentAmtLikelihoodT3 ==2).*ContractCorrect;
commitmentAmtLikelihood = sum(commitmentAmtLikelihoodT,2)./max(NumContractCorrect,1);

%%%%%%%% Don't want likelihood from commitment Amt to equal 0, otherwise log likelihood will evaluate to infinity. Keep small baseline probability here. 

commitmentAmtLikelihood = max(commitmentAmtLikelihood,0.1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Now genterate likelihood function for sequence of choices made by people in sample
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% including contract choice. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Here I will do this by matching entire sequence of choices predicted by candidate parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% To "Actual" choices made in simulated data. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Final Output of Likelihood Function 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%% People in control and control / discount only contribute through their Camp Attendence which is a binary variable
%%%%%%%%%%%%% NOTE: Go through these next five lines line by line, then tabulate Likelihood at each step to see what each type of person 
%%%%%%%%%%%%% is contributing to likelihood, given candidate set of parameters. 

T4A = (TreatmentPop==4).*(acceptsContractData(:,1)==1);
T4D = (TreatmentPop==4).*(acceptsContractData(:,1)==0);
T5A = (TreatmentPop==5).*(acceptsContractData(:,1)==1);
T5D = (TreatmentPop==5).*(acceptsContractData(:,1)==0);

%%%%%%%%%%%%%
%%%%%%%%%%%%%
%%%%%%%%%%%%%

Likelihood = zeros(Popsize,1);
Likelihood = Likelihood + Control.*attendsLikelihood + Discount.*attendsLikelihood;
Likelihood = Likelihood + T2.*acceptsContractLikelihood.*attendsLikelihood + T3.*acceptsContractLikelihood.*attendsLikelihood;
Likelihood = Likelihood + T4D.*acceptsContractLikelihood.*attendsLikelihood + T5D.*acceptsContractLikelihood.*attendsLikelihood; 
Likelihood = Likelihood + T4A.*acceptsContractLikelihood.*attendsLikelihood.*commitmentAmtLikelihood + T5A.*acceptsContractLikelihood.*attendsLikelihood.*commitmentAmtLikelihood; 
Likelihood = max(Likelihood, 0.01);

%%%%%%%%%%%%% Calculate Log-Likelihood, and take negative of this since algorithm in Matlab does minimizaiton, not maximization

output = -sum(log(Likelihood));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END OF FILE %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 